home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 5 / Amiga Tools 5.iso / grafik / 3d & render tools / irit / contrib / mrchcube / mc_main.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-07-16  |  13.2 KB  |  492 lines
  1. /******************************************************************************
  2. * An implementation of the marching cube algorithm (test file).              *
  3. *                                          *
  4. *                        Gershon Elber, Dec 1992.      *
  5. ******************************************************************************/
  6.  
  7. #include <stdio.h>
  8. #include "irit_sm.h"
  9. #include "mrchcube.h"
  10. #include "iritprsr.h"
  11. #include "allocate.h"
  12.  
  13. #define INPUT_ASCII    1
  14. #define INPUT_INTEGER    2
  15. #define INPUT_LONG    3
  16. #define INPUT_BYTE    4
  17. #define INPUT_FLOAT    5
  18. #define INPUT_DOUBLE    6
  19.  
  20.  
  21. static RealType
  22.     CubeWidth = 0.1;
  23. static int
  24.     DumpBinary = FALSE,
  25.     DumpOneLayer = 0,
  26.     DumpPolySkip = 1,
  27.     SkipInputData = 1,
  28.     DumpNormals = 1,
  29.     DataWidth = 100,
  30.     DataDepth = 100,
  31.     DataHeight = 100,
  32.     InputFormat = INPUT_ASCII;
  33.  
  34. static RealType GetOneScalar(void);
  35. static MCCubeCornerScalarStruct *GetCube(CagdBType FirstTime);
  36. static void EstimateGradient(MCCubeCornerScalarStruct *CCS);
  37. static void DumpPolys(MCPolygonStruct *Polys);
  38. static void DumpPolysBinary(MCPolygonStruct *Polys);
  39.  
  40. /******************************************************************************
  41. * Test routine of Marching Cubes.                          *
  42. * -t Threshold : the level (real number) to compute iso surface at.          *
  43. * -d Width Depth Height : dimension of sampled data.                  *
  44. * -f [ascii | int | long | byte | float | double] : single datum format.      *
  45. * -s n : dump out only every n polygons. 1 means dump all polygons.           *
  46. * -S n : skip n rows/cols/planes in input data. 1 means skip nothing.         *
  47. * -o n : process only layer number n.                          *
  48. * -n : do not do normals.                              *
  49. ******************************************************************************/
  50. void main(int argc, char **argv)
  51. {
  52.     MCCubeCornerScalarStruct *CCS;
  53.     RealType Threshold = 10.5;
  54.  
  55.     while (argc > 1) {
  56.     if (argv[1][0] == '-') {
  57.         switch (argv[1][1]) {
  58.         case 't':
  59.             sscanf(argv[2], "%lf", &Threshold);
  60.             fprintf(stderr, "Threshold set at %lf\n", Threshold);
  61.             argc -= 2;
  62.             argv += 2;
  63.             break;
  64.         case 'd':
  65.             sscanf(argv[2], "%d", &DataWidth);
  66.             sscanf(argv[3], "%d", &DataDepth);
  67.             sscanf(argv[4], "%d", &DataHeight);
  68.             fprintf(stderr, "Dimensions %dx%dx%d\n",
  69.                 DataWidth, DataDepth, DataHeight);
  70.             argc -= 4;
  71.             argv += 4;
  72.             break;
  73.         case 'f':
  74.             if (strncmp(argv[2], "ascii", 3) == 0)
  75.             InputFormat = INPUT_ASCII;
  76.             else if (strncmp(argv[2], "int", 3) == 0)
  77.             InputFormat = INPUT_INTEGER;
  78.             else if (strncmp(argv[2], "long", 3) == 0)
  79.             InputFormat = INPUT_LONG;
  80.             else if (strncmp(argv[2], "byte", 3) == 0)
  81.             InputFormat = INPUT_BYTE;
  82.             else if (strncmp(argv[2], "float", 3) == 0)
  83.             InputFormat = INPUT_FLOAT;
  84.             else if (strncmp(argv[2], "double", 3) == 0)
  85.             InputFormat = INPUT_DOUBLE;
  86.             else
  87.             {
  88.             fprintf(stderr, "Unknown input format \"%s\"\n",
  89.                 argv[2]);
  90.             exit(1);
  91.             }
  92.             fprintf(stderr, "Input format is set to \"%s\"\n",
  93.                 argv[2]);
  94.             argc -= 2;
  95.             argv += 2;
  96.             break;
  97.         case 's':
  98.             sscanf(argv[2], "%d", &DumpPolySkip);
  99.             fprintf(stderr, "DumpPolySkip is %d\n", DumpPolySkip);
  100.             argc -= 2;
  101.             argv += 2;
  102.             break;
  103.         case 'S':
  104.             sscanf(argv[2], "%d", &SkipInputData);
  105.             fprintf(stderr, "SkipInputData is %d\n", SkipInputData);
  106.             CubeWidth *= SkipInputData;
  107.             argc -= 2;
  108.             argv += 2;
  109.             break;
  110.         case 'n':
  111.             DumpNormals = 0;
  112.             fprintf(stderr, "Normal are not dumped\n");
  113.             argc--;
  114.             argv++;
  115.             break;
  116.         case 'b':
  117.             DumpBinary = TRUE;
  118.             fprintf(stderr, "Binary data is dumped out\n");
  119.             argc--;
  120.             argv++;
  121.             break;
  122.         case 'o':
  123.             sscanf(argv[2], "%d", &DumpOneLayer);
  124.             fprintf(stderr, "DumpOneLayer is %d\n", DumpOneLayer);
  125.             argc -= 2;
  126.             argv += 2;
  127.             break;
  128.         default:
  129.             fprintf(stderr, "Unknown option %c\n", argv[1][1]);
  130.             exit(1);
  131.         }
  132.     }
  133.     else {
  134.         fprintf(stderr, "Unknown command line %s\n", argv[1]);
  135.         exit(1);
  136.     }
  137.     }
  138.  
  139.     GetCube(TRUE);
  140.  
  141.     if (DumpBinary) {
  142.     MCPolygonStruct
  143.         *AllPolys = NULL;
  144.  
  145.     while (CCS = GetCube(FALSE)) {
  146.         MCPolygonStruct *PolyTmp,
  147.             *Polys = MCThresholdCube(CCS, Threshold);
  148.  
  149.         if (Polys != NULL) {
  150.         for (PolyTmp = Polys;
  151.              PolyTmp -> Pnext != NULL;
  152.              PolyTmp = PolyTmp -> Pnext);
  153.         PolyTmp -> Pnext = AllPolys;
  154.         AllPolys = Polys;
  155.         }
  156.     }
  157.  
  158.     DumpPolysBinary(AllPolys);
  159.     }
  160.     else {
  161.     printf("[OBJECT MC\n");
  162.     while (CCS = GetCube(FALSE)) {
  163.         MCPolygonStruct
  164.             *Polys = MCThresholdCube(CCS, Threshold);
  165.  
  166.         DumpPolys(Polys);
  167.     }
  168.     printf("]\n");
  169.     }
  170.  
  171.     exit(0);
  172. }
  173.  
  174. /******************************************************************************
  175. * Read stdin and returns one real at a time.                      *
  176. ******************************************************************************/
  177. static RealType GetOneScalar(void)
  178. {
  179.     short i;
  180.     long l;
  181.     float f;
  182.     double d;
  183.     RealType r;
  184.  
  185.     switch (InputFormat) {
  186.     case INPUT_ASCII:
  187.         if (scanf("%lf", &r) != 1)
  188.         return INFINITY;
  189.         break;
  190.     case INPUT_INTEGER:
  191.         i = getchar();
  192.         r = getchar() * 256 + i;
  193.         break;
  194.     case INPUT_LONG:
  195.         if (read(0, &l, 4) != 4)
  196.         return INFINITY;
  197.         r = l;
  198.         break;
  199.     case INPUT_BYTE:
  200.         if ((r = getchar()) == EOF)
  201.         return INFINITY;
  202.         break;
  203.     case INPUT_FLOAT:
  204.         if (read(0, &f, 4) != 4)
  205.         return INFINITY;
  206.         r = f;
  207.         break;
  208.     case INPUT_DOUBLE:
  209.         if (read(0, &d, 8) != 8)
  210.         return INFINITY;
  211.         r = d;
  212.         break;
  213.     default:
  214.         fprintf(stderr, "Input format requested not supported.\n");
  215.         exit(1);
  216.         break;
  217.     }
  218.  
  219.     return r;
  220. }
  221.  
  222. /******************************************************************************
  223. * Read stdin and returns one cube at a time.                      *
  224. ******************************************************************************/
  225. static MCCubeCornerScalarStruct *GetCube(CagdBType FirstTime)
  226. {
  227.     static MCCubeCornerScalarStruct
  228.     CCS;
  229.     static int
  230.     ProcessedOneLayer = FALSE,
  231.     LayerCountX = -1,
  232.     LayerCountY = 0,
  233.     LayerNumber = -1;
  234.     static RealType
  235.     *LayerOne = NULL,
  236.     *LayerTwo = NULL;
  237.     int i, j;
  238.     RealType *p;
  239.  
  240.     if (FirstTime) {
  241.     /* Initialize the CCS constant data */
  242.     CCS.CubeDim[0] = CubeWidth;
  243.     CCS.CubeDim[1] = CubeWidth;
  244.     CCS.CubeDim[2] = CubeWidth;
  245.     LayerNumber = -SkipInputData;
  246.     LayerCountX = -1;
  247.     LayerCountY = 0;
  248.  
  249.     if (LayerOne == NULL) {
  250.         LayerOne = (RealType *) IritMalloc(sizeof(RealType) * DataWidth *
  251.                            DataDepth);
  252.         LayerTwo = (RealType *) IritMalloc(sizeof(RealType) * DataWidth *
  253.                            DataDepth);
  254.     }
  255.  
  256.     for (p = LayerTwo, i = 0; i < DataWidth * DataDepth; i++)
  257.         if ((*p++ = GetOneScalar()) == INFINITY)
  258.         return NULL;
  259.     return NULL;
  260.     }
  261.  
  262.     if (LayerCountX == -1) {            /* Read the next layer. */
  263.     do {
  264.         LayerNumber += SkipInputData;
  265.         if ((DumpOneLayer > 0 && ProcessedOneLayer) ||
  266.         LayerNumber >= DataHeight - 1) {
  267.         IritFree((VoidPtr) LayerOne);
  268.         IritFree((VoidPtr) LayerTwo);
  269.         return NULL;
  270.         }
  271.  
  272.         p = LayerOne;
  273.         LayerOne = LayerTwo;
  274.         LayerTwo = p;
  275.  
  276.         for (j = 0; j < SkipInputData; j++) {
  277.         for (p = LayerTwo, i = 0; i < DataWidth * DataDepth; i++)
  278.             if ((*p++ = GetOneScalar()) == INFINITY)
  279.                 return NULL;
  280.         }
  281.  
  282.         LayerCountX = 0;
  283.         LayerCountY = 0;
  284.  
  285.         fprintf(stderr, "Doing layer %d\n", LayerNumber);
  286.     }
  287.     while (DumpOneLayer > LayerNumber);
  288.  
  289.     ProcessedOneLayer = TRUE;
  290.     }
  291.     
  292.  
  293.     CCS.Vrtx0Lctn[0] = LayerCountX * CubeWidth / SkipInputData;
  294.     CCS.Vrtx0Lctn[1] = LayerCountY * CubeWidth / SkipInputData;
  295.     CCS.Vrtx0Lctn[2] = LayerNumber * CubeWidth / SkipInputData;
  296.     CCS.Corners[0] = LayerOne[LayerCountY * DataWidth + LayerCountX];
  297.     CCS.Corners[1] = LayerOne[LayerCountY * DataWidth +
  298.                   LayerCountX + SkipInputData];
  299.     CCS.Corners[2] = LayerOne[(LayerCountY + SkipInputData) * DataWidth +
  300.                   LayerCountX + SkipInputData];
  301.     CCS.Corners[3] = LayerOne[(LayerCountY + SkipInputData) * DataWidth +
  302.                   LayerCountX];
  303.     CCS.Corners[4] = LayerTwo[LayerCountY * DataWidth + LayerCountX];
  304.     CCS.Corners[5] = LayerTwo[LayerCountY * DataWidth +
  305.                   LayerCountX + SkipInputData];
  306.     CCS.Corners[6] = LayerTwo[(LayerCountY + SkipInputData) * DataWidth +
  307.                   LayerCountX + SkipInputData];
  308.     CCS.Corners[7] = LayerTwo[(LayerCountY + SkipInputData) * DataWidth +
  309.                   LayerCountX];
  310.  
  311.     EstimateGradient(&CCS);
  312.  
  313.     LayerCountX += SkipInputData;
  314.     if (LayerCountX >= DataWidth - SkipInputData) {
  315.     LayerCountY += SkipInputData;
  316.     LayerCountX = 0;
  317.     if (LayerCountY >= DataDepth - SkipInputData)
  318.         LayerCountX = -1; /* No more data */
  319.     }
  320.  
  321.     return &CCS;
  322. }
  323.  
  324. /******************************************************************************
  325. * Estimate the gradient at the eight vertices, by first order difference.     *
  326. * Note gradient is estimated from this cube's eight vertices only.          *
  327. ******************************************************************************/
  328. static void EstimateGradient(MCCubeCornerScalarStruct *CCS)
  329. {
  330.     CCS -> GradientX[0] =
  331.     CCS -> GradientX[1] =
  332.         CCS -> Corners[1] - CCS -> Corners[0];
  333.     CCS -> GradientX[2] =
  334.     CCS -> GradientX[3] =
  335.         CCS -> Corners[2] - CCS -> Corners[3];
  336.     CCS -> GradientX[4] =
  337.     CCS -> GradientX[5] =
  338.         CCS -> Corners[5] - CCS -> Corners[4];
  339.     CCS -> GradientX[6] =
  340.     CCS -> GradientX[7] =
  341.         CCS -> Corners[6] - CCS -> Corners[7];
  342.     
  343.     CCS -> GradientY[0] =
  344.     CCS -> GradientY[3] =
  345.         CCS -> Corners[3] - CCS -> Corners[0];
  346.     CCS -> GradientY[1] =
  347.     CCS -> GradientY[2] =
  348.         CCS -> Corners[2] - CCS -> Corners[1];
  349.     CCS -> GradientY[4] =
  350.     CCS -> GradientY[7] =
  351.         CCS -> Corners[7] - CCS -> Corners[4];
  352.     CCS -> GradientY[5] =
  353.     CCS -> GradientY[6] =
  354.         CCS -> Corners[6] - CCS -> Corners[5];
  355.     
  356.     CCS -> GradientZ[0] =
  357.     CCS -> GradientZ[4] =
  358.         CCS -> Corners[4] - CCS -> Corners[0];
  359.     CCS -> GradientZ[1] =
  360.     CCS -> GradientZ[5] =
  361.         CCS -> Corners[5] - CCS -> Corners[1];
  362.     CCS -> GradientZ[2] =
  363.     CCS -> GradientZ[6] =
  364.         CCS -> Corners[6] - CCS -> Corners[2];
  365.     CCS -> GradientZ[3] =
  366.     CCS -> GradientZ[7] =
  367.         CCS -> Corners[7] - CCS -> Corners[3];
  368.  
  369.     CCS -> HasGradient = TRUE;
  370. }
  371.  
  372. /******************************************************************************
  373. * Dumps the polygons to stdout (Only triangles).                  *
  374. ******************************************************************************/
  375. static void DumpPolys(MCPolygonStruct *Polys)
  376. {
  377.     static int
  378.     DumpPolySkipLcl = 1;
  379.     char *VrtxFrmt;
  380.  
  381.     if (DumpNormals)
  382.     VrtxFrmt = "\t[[NORMAL %9.6lg %9.6lg %9.6lg] %9.6lg %9.6lg %9.6lg]\n";
  383.     else
  384.     VrtxFrmt = "\t[%9.6lg %9.6lg %9.6lg]\n";
  385.  
  386.     while (Polys) {
  387.     int i;
  388.     MCPolygonStruct
  389.         *Poly = Polys;
  390.  
  391.     Polys = Polys -> Pnext;
  392.  
  393.     /* Skip polygons in the output (so we can display them - they      */
  394.     /* are just too many of them.                       */
  395.     if (--DumpPolySkipLcl <= 0) {
  396.         DumpPolySkipLcl = DumpPolySkip;
  397.  
  398.         for (i = 2; i < Poly -> NumOfVertices - 1; i++) {
  399.         printf("    [POLYGON 3\n");
  400.         if (DumpNormals) {
  401.             printf(VrtxFrmt,
  402.                Poly -> N[0][0], Poly -> N[0][1], Poly -> N[0][2],
  403.                Poly -> V[0][0], Poly -> V[0][1], Poly -> V[0][2]);
  404.             printf(VrtxFrmt,
  405.                Poly -> N[i-1][0], Poly -> N[i-1][1], Poly -> N[i-1][2],
  406.                Poly -> V[i-1][0], Poly -> V[i-1][1], Poly -> V[i-1][2]);
  407.             printf(VrtxFrmt,
  408.                Poly -> N[i][0], Poly -> N[i][1], Poly -> N[i][2],
  409.                Poly -> V[i][0], Poly -> V[i][1], Poly -> V[i][2]);
  410.         }
  411.         else {
  412.             printf(VrtxFrmt,
  413.                Poly -> V[0][0], Poly -> V[0][1], Poly -> V[0][2]);
  414.             printf(VrtxFrmt,
  415.                Poly -> V[i-1][0], Poly -> V[i-1][1], Poly -> V[i-1][2]);
  416.             printf(VrtxFrmt,
  417.                Poly -> V[i][0], Poly -> V[i][1], Poly -> V[i][2]);
  418.         }
  419.         printf("    ]\n");
  420.         }
  421.     }
  422.  
  423.     IritFree((VoidPtr) Poly);
  424.     }
  425. }
  426.  
  427. /******************************************************************************
  428. * Dumps triangles to stdout as binary file.                      *
  429. ******************************************************************************/
  430. static void DumpPolysBinary(MCPolygonStruct *Polys)
  431. {
  432.     static int
  433.     DumpPolySkipLcl = 1;
  434.     int i, j, Handler;
  435.     IPPolygonStruct *Pl,
  436.     *PlHead = NULL;
  437.     IPObjectStruct *PObj;
  438.  
  439.     while (Polys) {
  440.     IPVertexStruct *V1, *V2, *V3;
  441.     int i;
  442.     MCPolygonStruct
  443.         *Poly = Polys;
  444.  
  445.     Polys = Polys -> Pnext;
  446.  
  447.     /* Skip polygons in the output (so we can display them - they      */
  448.     /* are just too many of them.                       */
  449.     if (--DumpPolySkipLcl <= 0) {
  450.         DumpPolySkipLcl = DumpPolySkip;
  451.  
  452.         for (i = 2; i < Poly -> NumOfVertices - 1; i++) {
  453.         V3 = IPAllocVertex(0, 0, NULL, NULL);
  454.         V2 = IPAllocVertex(0, 0, NULL, V3);
  455.         V1 = IPAllocVertex(0, 0, NULL, V2);
  456.  
  457.         Pl = IPAllocPolygon(0, 0, V1, PlHead);
  458.         PlHead = Pl;
  459.  
  460.         for (j = 0; j < 3; j++) {
  461.             V1 -> Coord[j] = Poly -> V[0][j];
  462.             V2 -> Coord[j] = Poly -> V[i - 1][j];
  463.             V3 -> Coord[j] = Poly -> V[i][j];
  464.         }
  465.  
  466.         if (DumpNormals) {
  467.             for (j = 0; j < 3; j++) {
  468.             V1 -> Normal[j] = Poly -> N[0][j];
  469.             V2 -> Normal[j] = Poly -> N[i - 1][j];
  470.             V3 -> Normal[j] = Poly -> N[i][j];
  471.             }
  472.             IP_SET_NORMAL_VRTX(V1);
  473.             IP_SET_NORMAL_VRTX(V2);
  474.             IP_SET_NORMAL_VRTX(V3);
  475.         }
  476.         }
  477.     }
  478.  
  479.     IritFree((VoidPtr) Poly);
  480.     }
  481.  
  482.     PObj = IPAllocObject("MCH", IP_OBJ_POLY, NULL);
  483.     PObj -> U.Pl = PlHead;
  484.  
  485.     Handler = IritPrsrOpenStreamFromFile(stdout, FALSE, TRUE, FALSE);
  486.     IritPrsrPutObjectToHandler(Handler, PObj);
  487.     IritPrsrCloseStream(Handler, TRUE);
  488.  
  489.     IPFreeObject(PObj);
  490. }
  491.  
  492.